Cours_Geneve_6

JE prépare ma session de travail

setwd("~/GitHub/Cours_2020_UniGE/Cours_Geneve_6")

Je charge les packages utiles

if(!require("stylo")){
  install.packages("stylo")
  library(stylo)
}
if(!require("cluster")){
  install.packages("cluster")
  library(cluster)
}
if(!require("FactoMineR")){
  install.packages("FactoMineR")
  library(FactoMineR)
}
if(!require("factoextra")){
  install.packages("factoextra")
  library(factoextra)
}
if(!require("smacof")){
  install.packages("smacof")
  library(smacof)
}

Loading required package: smacof
Loading required package: plotrix
Registered S3 methods overwritten by 'car':
  method                          from
  influence.merMod                lme4
  cooks.distance.influence.merMod lme4
  dfbeta.influence.merMod         lme4
  dfbetas.influence.merMod        lme4

Attaching package: ‘smacof’

The following object is masked from ‘package:base’:

    transform

Ma première CAH depuis 0

Une classification ascendante hiérarchique.

# Je crée un tableau avec des données au hasard
tableau <- cbind(x=c(1,2,5,8,11),y=c(2,4,1,12,10))
colnames(tableau)<- c("Avoir", "Être")
rownames(tableau) <- c("Racine","Molière","Racine","Scarron","Corneille")
tableau

          Avoir Être
Racine        1    2
Molière       2    4
Racine        5    1
Scarron       8   12
Corneille    11   10

Calcul de distance, par défaut eucidienne. on peut utiliser le paramètre method=“manhattan” pour une autre distance

distance_tableau <- dist(tableau)
distance_tableau

             Racine   Molière    Racine   Scarron
Molière    2.236068                              
Racine     4.123106  4.242641                    
Scarron   12.206556 10.000000 11.401754          
Corneille 12.806248 10.816654 10.816654  3.605551

Classification/clustering. Par défaut complete linkage. On peut utiliser le paramètre method=“centroid” ou “ward.D”

hc <- hclust(distance_tableau)
plot(hc)

#On customise le graph
plot(hc, main="Ma première classification",  xlab="Euclidian distance, \n Complete-linkage clustering" , sub="" , hang =-1)

On peut faire soi-même un delta de Burrows. Pour cela je z-score mes données

#Je transforme mon tableau en data.frame
tableau_df<-as.data.frame(tableau)
#Je z-score mes deux colonnes avec la fonction scale()
tableau_z_a<-scale(tableau_df$Avoir)
tableau_z_b<-scale(tableau_df$Être)
#Je reforme un tableau, avec les nouvelles données
tableau_z<-cbind(tableau_z_a,tableau_z_b)
tableau_z

            [,1]       [,2]
[1,] -1.05786348 -0.7724598
[2,] -0.81743996 -0.3659020
[3,] -0.09616941 -0.9757388
[4,]  0.62510115  1.2603292
[5,]  1.34637170  0.8537714

Et je produis un nouveau graph

hc <- hclust(distance_tableau, method="ward.D")
plot(hc, main="Ma première classificatoin",  xlab="Burrows Delta, \n Ward's method" , sub="" , hang =-1)

Un premier cas

Je charge mon corpus

BOYER_AMOURSJUPITERSEMELE_1666<-paste(scan("corpus/BOYER_AMOURSJUPITERSEMELE_1666.txt", what="character", sep=""),collapse=" ")

Read 17734 items

BOYER_ARISTODEME_1648<-paste(scan("corpus/BOYER_ARISTODEME_1648.txt", what="character", sep=""),collapse=" ")

Read 11106 items

CORNEILLEP_ANDROMEDE_1651<-paste(scan("corpus/CORNEILLEP_ANDROMEDE_1651.txt", what="character", sep=""),collapse=" ")

Read 16121 items

CORNEILLEP_MEDEE_1682<-paste(scan("corpus/CORNEILLEP_MEDEE_1682.txt", what="character", sep=""),collapse=" ")

Read 13705 items

DURYER_DYNAMIS_1653<-paste(scan("corpus/DURYER_DYNAMIS_1653.txt", what="character", sep=""),collapse=" ")

Read 16793 items

DURYER_CLITOPHON<-paste(scan("corpus/DURYER_CLITOPHON_1632.txt", what="character", sep=""),collapse=" ")

Read 15265 items

MOLIERE_AMPHITRYON_1668<-paste(scan("corpus/MOLIERE_AMPHITRYON_1668.txt", what="character", sep=""),collapse=" ")

Read 14152 items

MOLIERE_MISANTHROPE_1667<-paste(scan("corpus/MOLIERE_MISANTHROPE_1667.txt", what="character", sep=""),collapse=" ")

Read 16074 items

RACINE_PHEDRE_1697<-paste(scan("corpus/RACINE_PHEDRE_1697.txt", what="character", sep=""),collapse=" ")

Read 13434 items

RACINE_ESTHER_1689<-paste(scan("corpus/RACINE_ESTHER_1689.txt", what="character", sep=""),collapse=" ")

Read 10460 items

SCARRON_DOMJAPHETDARMENIE_1653<-paste(scan("corpus/SCARRON_DOMJAPHETDARMENIE_1653.txt", what="character", sep=""),collapse=" ")

Read 12945 items

SCARRON_ECOLIERDESALAMANQUE_1655<-paste(scan("corpus/SCARRON_ECOLIERDESALAMANQUE_1655.txt", what="character", sep=""),collapse=" ")

Read 14133 items

ROTROU_HERCULEMOURANT_1636<-paste(scan("corpus/ROTROU_HERCULEMOURANT_1636.txt", what="character", sep=""),collapse=" ")

Read 12561 items

ROTROU_DOMBERNARDDECABRERE_1647<-paste(scan("corpus/ROTROU_DOMBERNARDDECABRERE_1647.txt", what="character", sep=""),collapse=" ")

Read 14546 items

SCUDERY_MORTDECESAR_1637<-paste(scan("corpus/SCUDERY_MORTDECESAR_1637.txt", what="character", sep=""),collapse=" ")

Read 11918 items

SCUDERY_ORANTE_1636<-paste(scan("corpus/SCUDERY_ORANTE_1636.txt", what="character", sep=""),collapse=" ")

Read 12408 items

#Je peux regarder ici un exemple (enlever le # au début de la ligne)
#CORNEILLEP_ANDROMEDE_1651

J’en fais un tableau

#Je crée une liste de mes textes
my.corpus.raw = list(BOYER_AMOURSJUPITERSEMELE_1666,
                     BOYER_ARISTODEME_1648,
                     DURYER_DYNAMIS_1653,
                     DURYER_CLITOPHON,
                     MOLIERE_AMPHITRYON_1668,
                     MOLIERE_MISANTHROPE_1667,
                     SCARRON_DOMJAPHETDARMENIE_1653,
                     SCARRON_ECOLIERDESALAMANQUE_1655)
#Je tokenise mon texte (espace, apostrophe, retrait de la ponctuation)
my.corpus.clean = lapply(my.corpus.raw, txt.to.words)
#Je compte la fréquence des tokens
complete.word.list = make.frequency.list(my.corpus.clean)
#Je transforme le résultat en table, où les résultats de chaque mot/texte sont alignés
table.of.frequencies=make.table.of.frequencies(my.corpus.clean, complete.word.list, relative = F)

processing  8  text samples
       
combining frequencies into a table...

#Je donne un nom à mes colonnes
row.names(table.of.frequencies)=c("boyer_agamemnon",
                                  "boyer_amoursJupiter",
                                  "duryer_dynamis",
                                  "duryer_clitophon",
                                  "moliere_amphytrion",
                                  "moliere_misanthrope",
                                  "scarron_domJaphet",
                                  "scarron_ecolierSalamanque")
#Je sauvegarde une copie
write.csv(table.of.frequencies, file = "table.of.frequencies.csv",row.names=TRUE)
#je convertis mes données en dataframe
table.of.frequencies = as.data.frame(read.csv(file="table.of.frequencies.csv", sep = ",", header = TRUE, row.names=1, quote = '\"'))
#J'affiche le résultat
#View(table.of.frequencies)

Ce résultat ne m’arrange pas: je vais avoir besoin d’avoir les occurrences en rang, et les textes en colonne, et donc d’inverser mon tableau

#j'utilise la fonction transpose(?)
table.of.frequencies<-t(table.of.frequencies)
#View(table.of.frequencies)

J’observe la distribution de mon corpus.

summary(table.of.frequencies)

 boyer_agamemnon   boyer_amoursJupiter duryer_dynamis    duryer_clitophon moliere_amphytrion
 Min.   :  0.000   Min.   :  0.000     Min.   :  0.000   Min.   :  0.00   Min.   :  0.000   
 1st Qu.:  0.000   1st Qu.:  0.000     1st Qu.:  0.000   1st Qu.:  0.00   1st Qu.:  0.000   
 Median :  0.000   Median :  0.000     Median :  0.000   Median :  0.00   Median :  0.000   
 Mean   :  2.208   Mean   :  1.425     Mean   :  2.127   Mean   :  1.92   Mean   :  1.783   
 3rd Qu.:  1.000   3rd Qu.:  0.000     3rd Qu.:  1.000   3rd Qu.:  1.00   3rd Qu.:  1.000   
 Max.   :599.000   Max.   :321.000     Max.   :589.000   Max.   :537.00   Max.   :478.000   
 moliere_misanthrope scarron_domJaphet scarron_ecolierSalamanque
 Min.   :  0.000     Min.   :  0.000   Min.   :  0.000          
 1st Qu.:  0.000     1st Qu.:  0.000   1st Qu.:  0.000          
 Median :  0.000     Median :  0.000   Median :  0.000          
 Mean   :  2.018     Mean   :  1.686   Mean   :  1.828          
 3rd Qu.:  1.000     3rd Qu.:  1.000   3rd Qu.:  1.000          
 Max.   :567.000     Max.   :425.000   Max.   :423.000          

Je peux synthétiser ce résultat avec une boîte à moustache des fréquences (absolues) des mots du corpus.

100% center

Source: stat4decision

maBoite <- boxplot(colSums(table.of.frequencies), main = "Distribution du nombre de mots par de texte", ylab="Nombre de mots (fréq. absolues)", sub="Corpus: distribution")

Je contrôle l’hétérogénéité de mon corpus

#Je dois regrouper les données par auteur
BOYER = colSums(table.of.frequencies[,grepl('boyer_', colnames(table.of.frequencies))])
CORNEILLE = colSums(table.of.frequencies[,grepl('corneille_', colnames(table.of.frequencies))])
DURYER = colSums(table.of.frequencies[,grepl('duryer_', colnames(table.of.frequencies))])
MOLIERE = colSums(table.of.frequencies[,grepl('moliere_', colnames(table.of.frequencies))])
RACINE = colSums(table.of.frequencies[,grepl('racine_', colnames(table.of.frequencies))])
SCARRON = colSums(table.of.frequencies[,grepl('scarron_', colnames(table.of.frequencies))])
ROTROU = colSums(table.of.frequencies[,grepl('rotrou_', colnames(table.of.frequencies))])
#Et le graphique
boxplot(list(BOYER,MOLIERE, CORNEILLE,DURYER,RACINE,SCARRON, ROTROU), names=c('BOYER','MOLIERE', 'CORNEILLE','DURYER','RACINE','SCARRON','ROTROU'), main="Longueur des textes", sub="Corpus: distribution",ylab="Nombre de mots (fréquences absolues)")

Je transforme cette table des fréquences absolues en une table des fréquences relatives

#Je fais une copie de ma table de fréquences
freqs_rel = table.of.frequencies
#Dans cette copie, pour chacun des mots de chaque colonne, je divise le chiffre trouvé par la somme de la colonne
for(i in 1:ncol(freqs_rel)){
        freqs_rel[,i] = freqs_rel[,i]/sum(freqs_rel[,i])
}
head(freqs_rel)

     boyer_agamemnon boyer_amoursJupiter duryer_dynamis duryer_clitophon moliere_amphytrion
de        0.03211280          0.02666556     0.02459929       0.03310319         0.03173339
et        0.02728784          0.02475494     0.03278050       0.02404143         0.02701985
je        0.01613681          0.01918923     0.01775378       0.01658242         0.02243909
vous      0.02299898          0.01844160     0.01703028       0.00961657         0.01566753
que       0.01704820          0.01910616     0.02276269       0.02416471         0.02077939
le        0.01715542          0.01702941     0.02075913       0.01713722         0.01427339
     moliere_misanthrope scarron_domJaphet scarron_ecolierSalamanque
de            0.03162403        0.02984131                0.02738573
et            0.03150669        0.02260918                0.02732099
je            0.02434874        0.02569864                0.02596141
vous          0.03326684        0.02155596                0.02537874
que           0.02018306        0.01530684                0.02026415
le            0.01513729        0.01987080                0.01573223

Je peux (et même je dois) sélectionner dans cette liste les n plus fréquents

#Je prends les premiers rangs. Changez ce chiffre pour changer les résultats par la suite
freqs_rel_mfw = freqs_rel[1:100,]

Analyse

La CAH

Regardons à quoi ressemble notre CAH

CAH = agnes(as.dist(dist.wurzburg(t(freqs_rel_mfw))), method="ward")
CAH_orig_aggloCoeff <- CAH$ac
plot(CAH)

Nettoyons tout cela…

plot(CAH,
     main="Cluster analysis",
     xlab=paste("Cosine Delta\n Agglomerative coefficient =", CAH_orig_aggloCoeff),
     hang = -0.1)

Faisons un peu mieux

#on va colorier les vecteurs en fonction de l'auteur (en pratique, la chaîne de caractère avant l'_underscore_)
labels_color = vector(length = length(CAH$order.lab))
labels_color[grep("boyer", CAH$order.lab)] = "darkblue"
labels_color[grep("duryer", CAH$order.lab)] = "darkgrey"
labels_color[grep("scarron", CAH$order.lab)] = "darkred"
labels_color[grep("moliere", CAH$order.lab)] = "darkgreen"
#je reprends mon graph
plot(CAH)

#Je fais un peu de coloriage
factoextra::fviz_dend(CAH,
                      ylab=paste("Cosine Delta\n Agglomerative coefficient =", round(CAH_orig_aggloCoeff, digit=3)),
                      k=2, # je cherche à dégager deux clusters
                      rect = TRUE,
                      k_colors = rep("black",2), #je garde la couleur des traits en noir
                      labels_track_height = 0.7, # Je gère la taille de mon rectangle
                      label_cols = labels_color, #j'applique mes couleurs
                      horiz=T) # je mets mon graph à l'horizontal, parce que c'est plus pratique

C’est bien mieux, mais pas encore parfait: le code est est trop long, pas très propre. Cachons tout cela dans un coin. Pour cela, compliquons-nous le travail un peu plus en créant une fonction je que conserve dans un autre fichier appelé functions.R.

Je vais donc charger cette fonction avant de l’utiliser.

source("functions.R")

Quatre paramètres ont été prévus * Le nom des données * La légende à afficher * Le coeffichient d’agglomération (récupéré au passage, un peu plus haut) * une variable pour le paramètre labels_track_height

customPlot(CAH,"Original texts, 100 MFW, Culled @ 0%, Distance: wurzburg",CAH_orig_aggloCoeff,0.7)

Je peux donc reproduire ma CAH simplement avec d’autres configurations. Tentons désormais avec plus de MFW

#on passe de 100 à 1000
freqs_rel_mfw = freqs_rel[1:1000,]
#je recalcule ma CAH
CAH = agnes(as.dist(dist.wurzburg(t(freqs_rel_mfw))), method="ward")
CAH_orig_aggloCoeff <- CAH$ac
#Et désormais, je n'ai plus qu'une simple fonction "maison" pour me produire ma visualisation
customPlot(CAH,"Original texts, 1000 MFW, Culled @ 0%, Distance: wurzburg",CAH_orig_aggloCoeff,0.5)

Tentons avec une autre distance: Ruzicka measure ou Minmax (cf.Koppel, M. and Winter, Y. (2014). Determining if two documents are written by the same author. Journal of the Association for Information Science and Technology, 65: 178–87) dont la robustesse a été démontrée par Kestmont et al. ( Kestemont, M.Stover, J.Koppel, M.Karsdorp, K.Daelemans, W., Authorship Verification with the Ruzicka Metric, DH 2016, Jul 2016, Krakow (Poland)).

100% center

(tf pour term frequency)

#je recalcule ma CAH
CAH = agnes(as.dist(dist.minmax(t(freqs_rel_mfw))), method="ward")
CAH_orig_aggloCoeff <- CAH$ac
#Et désormais, je n'ai plus qu'une simple fonction "maison" pour me produire ma visualisation
customPlot(CAH,"Original texts, 1000 MFW, Culled @ 0%, Distance: Minmax",CAH_orig_aggloCoeff,-0.2)

Et encore une autre, en diminuant de nouveau les MFW

#on repasse de 1000 à 100
freqs_rel_mfw = freqs_rel[1:100,]
#je recalcule ma CAH
CAH = agnes(as.dist(dist.wurzburg(t(freqs_rel_mfw))), method="ward")
CAH_orig_aggloCoeff <- CAH$ac
customPlot(CAH,"Original texts, 100 MFW, Culled @ 0%, Distance: wurzburg",CAH_orig_aggloCoeff,0.7)

Under the hood: les classes

Tout cela est très bien. Mais qu’est-ce qui est caché derrière ces clusters? Nous pouvons le savoir. Prenons notre dernière CAH, et regardons le processus de construction des clusters:

CAH_2 = as.hclust(CAH)
plot(CAH_2$height, type="h", ylab="height")

Je remarque que j’ai plusieurs clusters. Tentons de voir ce qu’ils contiennent. Classons par 4 notre corpus:

#Je "coupe en deux" ma CAH, ou plutôt je demande une classification en deux des différents
classes = cutree(CAH, k = "4")
#J'ajoute mes classes aux fréquences que j'ai utilisées
Processed_classes = t(freqs_rel_mfw)
Processed_classes = cbind(as.data.frame(Processed_classes), as.factor(classes))
colnames(Processed_classes)[101] = "Classes"
#J'affiche le résultat
Processed_classes[101]

                          Classes
boyer_agamemnon                 1
boyer_amoursJupiter             1
duryer_dynamis                  2
duryer_clitophon                2
moliere_amphytrion              3
moliere_misanthrope             3
scarron_domJaphet               4
scarron_ecolierSalamanque       4

Diminuons à deux le nombre de cluster

#Je "coupe en deux" ma CAH, ou plutôt je demande une classification en deux des différents
classes = cutree(CAH, k = "2")
#J'ajoute mes classes aux fréquences que j'ai utilisées
Processed_classes = t(freqs_rel_mfw)
Processed_classes = cbind(as.data.frame(Processed_classes), as.factor(classes))
colnames(Processed_classes)[101] = "Classes"
#J'affiche le résultat
Processed_classes[101]

                          Classes
boyer_agamemnon                 1
boyer_amoursJupiter             1
duryer_dynamis                  1
duryer_clitophon                1
moliere_amphytrion              2
moliere_misanthrope             2
scarron_domJaphet               2
scarron_ecolierSalamanque       2

Je peux aussi regarder les valeurs associées à chacun des mots, et ceux associés à chaque cluster:

myClasses = catdes(Processed_classes, num.var = 101)
myClasses

Link between the cluster variable and the quantitative variables
================================================================
           Eta2      P-value
je    0.8855821 0.0004896757
amour 0.8040975 0.0025455982
yeux  0.7601916 0.0047643342
moi   0.7512318 0.0053410502
y     0.7445267 0.0058031317
pas   0.7375033 0.0063162127
être  0.7371581 0.0063422162
n     0.7356576 0.0064561179
si    0.7193183 0.0077904004
l     0.7138624 0.0082757444
ne    0.6959504 0.0100182228
faire 0.6770794 0.0121175393
dieux 0.6744536 0.0124323540
qui   0.6573420 0.0146275157
mes   0.6342385 0.0180111247
en    0.6330403 0.0182005359
veux  0.6181419 0.0206765720
m     0.6095886 0.0222026475
ses   0.5960834 0.0247753466
plus  0.5890029 0.0262070694
suis  0.5741345 0.0294079634
comme 0.5457179 0.0363077270
bien  0.5423877 0.0371876525
me    0.5098090 0.0466527660
point 0.5027205 0.0489296469

Description of each cluster by quantitative variables
=====================================================
$`1`
         v.test Mean in category Overall mean sd in category   Overall sd    p.value
amour  2.372484      0.004818126  0.003380804   0.0008299353 0.0016028752 0.01766891
yeux   2.306803      0.001808802  0.001422762   0.0001256128 0.0004427618 0.02106579
si     2.243931      0.009360696  0.007472486   0.0015781275 0.0022263300 0.02483681
l      2.235405      0.017151679  0.015865132   0.0009675451 0.0015227141 0.02539075
dieux  2.172826      0.003528194  0.001978623   0.0014225437 0.0018868420 0.02979338
qui    2.145086      0.009527600  0.008502894   0.0005062741 0.0012638733 0.03194597
mes    2.107052      0.003316553  0.002528891   0.0007548847 0.0009890388 0.03511304
ses    2.042690      0.002753466  0.001970091   0.0007580942 0.0010146507 0.04108310
plus   2.030522      0.005931663  0.005243258   0.0006308718 0.0008969846 0.04230349
suis  -2.004730      0.001103213  0.002013334   0.0005036370 0.0012011368 0.04499194
m     -2.065701      0.004543850  0.005795718   0.0008642927 0.0016033934 0.03885672
veux  -2.080143      0.001166075  0.001530247   0.0001085254 0.0004631936 0.03751245
en    -2.105061      0.010588157  0.012344905   0.0015522742 0.0022079724 0.03528599
faire -2.177052      0.001783577  0.002343016   0.0004299235 0.0006798804 0.02947669
ne    -2.207182      0.008047953  0.009833099   0.0010665234 0.0021398561 0.02730134
n     -2.269274      0.005291597  0.006619240   0.0008650303 0.0015479011 0.02325169
être  -2.271587      0.001702773  0.002472563   0.0003633443 0.0008965862 0.02311147
pas   -2.272119      0.004991906  0.006090712   0.0006925025 0.0012794955 0.02307934
y     -2.282912      0.001141360  0.001708866   0.0003815097 0.0006577039 0.02243556
moi   -2.293169      0.003940317  0.005481745   0.0003737930 0.0017784282 0.02183829
je    -2.489794      0.017415563  0.021013767   0.0011821199 0.0038235906 0.01278171

$`2`
         v.test Mean in category Overall mean sd in category   Overall sd    p.value
je     2.489794     0.0246119713  0.021013767   0.0013957565 0.0038235906 0.01278171
moi    2.293169     0.0070231732  0.005481745   0.0011974515 0.0017784282 0.02183829
y      2.282912     0.0022763722  0.001708866   0.0002747232 0.0006577039 0.02243556
pas    2.272119     0.0071895170  0.006090712   0.0006163696 0.0012794955 0.02307934
être   2.271587     0.0032423534  0.002472563   0.0005390367 0.0008965862 0.02311147
n      2.269274     0.0079468818  0.006619240   0.0007200349 0.0015479011 0.02325169
ne     2.207182     0.0116182445  0.009833099   0.0012833566 0.0021398561 0.02730134
faire  2.177052     0.0029024541  0.002343016   0.0003371910 0.0006798804 0.02947669
en     2.105061     0.0141016520  0.012344905   0.0010809285 0.0022079724 0.03528599
veux   2.080143     0.0018944192  0.001530247   0.0003899695 0.0004631936 0.03751245
m      2.065701     0.0070475865  0.005795718   0.0011226719 0.0016033934 0.03885672
suis   2.004730     0.0029234554  0.002013334   0.0009875056 0.0012011368 0.04499194
plus  -2.030522     0.0045548539  0.005243258   0.0005131879 0.0008969846 0.04230349
ses   -2.042690     0.0011867150  0.001970091   0.0005069224 0.0010146507 0.04108310
mes   -2.107052     0.0017412296  0.002528891   0.0003817371 0.0009890388 0.03511304
qui   -2.145086     0.0074781878  0.008502894   0.0009156383 0.0012638733 0.03194597
dieux -2.172826     0.0004290517  0.001978623   0.0005425611 0.0018868420 0.02979338
l     -2.235405     0.0145785853  0.015865132   0.0006251135 0.0015227141 0.02539075
si    -2.243931     0.0055842772  0.007472486   0.0005403116 0.0022263300 0.02483681
yeux  -2.306803     0.0010367229  0.001422762   0.0002797223 0.0004427618 0.02106579
amour -2.372484     0.0019434819  0.003380804   0.0005637698 0.0016028752 0.01766891

Under the hood: l’es classes’ACP

Je commence par faire une ACP.

ACP = PCA(t(freqs_rel_mfw))

Faisons un peu de ménage, ajoutons de l’information, et tentons de comprendre

#categories
get_categories = vector(length = length(colnames(freqs_rel_mfw)))
get_categories[grep("racine", colnames(freqs_rel_mfw))] = "Racine"
get_categories[grep("moliere", colnames(freqs_rel_mfw))] = "Moliere"
get_categories[grep("duryer", colnames(freqs_rel_mfw))] = "Duryer"
get_categories[grep("boyer", colnames(freqs_rel_mfw))] = "Boyer"
ACP = PCA(t(freqs_rel_mfw))

fviz_pca_ind(ACP, col.ind = get_categories, legend="none")

On voit nettement deux axes, qu doivent bien s’appuyer sur quelque chose. Serait-il possible de faire ressortir les mots qui sont à l’origine de ce partitionnement des données?

Continuons notre nettoyage pour y voir plus clair

#categories
get_categories = vector(length = length(colnames(freqs_rel_mfw)))
get_categories[grep("racine", colnames(freqs_rel_mfw))] = "Racine"
get_categories[grep("moliere", colnames(freqs_rel_mfw))] = "Moliere"
ACP = PCA(t(freqs_rel_mfw))

fviz_pca_ind(ACP, col.ind = get_categories, legend="none")

fviz_pca_var(ACP, col.var="contrib",geom.var = "text", select.var = list(contrib =20))+
scale_color_gradient2(
                      low="yellow", mid="orange", 
                      high="darkred", midpoint=1.6)+theme_bw()

Nous nous rappelons qu’il existe un pourcentage de réalité retenu par axe, que nous pouvons observer, et donc calculer. Dans ce cas précis, il y en a seulement trois axes, car nous n’avons mis que quatre de textes.

ACP$eig

       eigenvalue percentage of variance cumulative percentage of variance
comp 1  30.037372              30.037372                          30.03737
comp 2  21.041028              21.041028                          51.07840
comp 3  15.055366              15.055366                          66.13377
comp 4  11.101088              11.101088                          77.23485
comp 5   9.201767               9.201767                          86.43662
comp 6   7.811868               7.811868                          94.24849
comp 7   5.751511               5.751511                         100.00000

Un petit graph montre bien la lente perte de significativité des axes :

barplot(ACP$eig[,1], main="Percentage of variance", names.arg=1:nrow(ACP$eig))

Les autres algorithmes

t-SNE

Avec l’algorithme de Barnes-Hut,

library(Rtsne)
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 0.5, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

Attention! Si on relance la même commande, on obtient un résultat différent!

library(Rtsne)
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 0.5, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

#encore une fois
library(Rtsne)
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 0.5, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

#encore une fois
library(Rtsne)
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 0.5, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

On peut faire varier la perpléxité (c’est à dire s’autoriser une distorsion plus grande): on passe à 1.5

maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 1.5, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

#Encore une fois
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 1.5, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

#Encore une fois
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 1.5, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

#Encore une fois
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 1.5, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

On peut faire varier la perpléxité: on passe à 2

maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 2, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

#Encore une fois
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 2, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

#Encore une fois
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 2, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

#Encore une fois
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 2, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

MDS

Le Multidimensional scaling est une compression des données en deux dimensions. Il en existe trois types: * classique * métrique * non métrique

MDS classique

Le MDS classique se base sur le calcul de distance, dont il va essayer de préserver l’essentiel.

#Je crée mes données
fit = cmdscale(dist(t(freqs_rel_mfw), method = "manhattan"), eig=TRUE, k=2) # k est le nombre de dimensions souhaité
#Je dessine mon résultat
x = fit$points[,1]
y = fit$points[,2]
plot(x, y, xlab=paste("Coordinate 1 (GOF: ", round(fit$GOF[1] * 100, digits=2), "%)"), ylab=paste("Coordinate 2 (GOF: ", round(fit$GOF[2] * 100, digits=2), "%)"), main="PMD métrique")
text(x, y, labels = row.names(t(freqs_rel_mfw)), cex=.7) 

MDS Métrique

Le MDS métrique, dit aussi ordinal, ne s’intéresse pas à la mesure de distance, mais sa valeur en relation avec les autres paires

Rappelons que le stress permet de mesurer la distortion introduite dans le résultat.

#J'affiche le stress avec le paramètre verbose=T(RUE)
MDSmetrique = mds(dist(t(freqs_rel_mfw), method = "manhattan"), ndim=2, type="ordinal",verbose=T)

Iteration:    1  Stress (raw):   0.01510422  Difference:    0.06233185 
Iteration:    2  Stress (raw):   0.00887573  Difference:    0.00622848 
Iteration:    3  Stress (raw):   0.00656184  Difference:    0.00231389 
Iteration:    4  Stress (raw):   0.00541200  Difference:    0.00114984 
Iteration:    5  Stress (raw):   0.00475869  Difference:    0.00065331 
Iteration:    6  Stress (raw):   0.00439233  Difference:    0.00036636 
Iteration:    7  Stress (raw):   0.00417426  Difference:    0.00021807 
Iteration:    8  Stress (raw):   0.00403217  Difference:    0.00014209 
Iteration:    9  Stress (raw):   0.00393218  Difference:    0.00009999 
Iteration:   10  Stress (raw):   0.00385734  Difference:    0.00007485 
Iteration:   11  Stress (raw):   0.00379958  Difference:    0.00005776 
Iteration:   12  Stress (raw):   0.00375696  Difference:    0.00004262 
Iteration:   13  Stress (raw):   0.00372424  Difference:    0.00003271 
Iteration:   14  Stress (raw):   0.00369790  Difference:    0.00002634 
Iteration:   15  Stress (raw):   0.00367617  Difference:    0.00002173 
Iteration:   16  Stress (raw):   0.00365798  Difference:    0.00001819 
Iteration:   17  Stress (raw):   0.00364262  Difference:    0.00001537 
Iteration:   18  Stress (raw):   0.00362954  Difference:    0.00001308 
Iteration:   19  Stress (raw):   0.00361836  Difference:    0.00001119 
Iteration:   20  Stress (raw):   0.00360875  Difference:    0.00000961 
Iteration:   21  Stress (raw):   0.00360047  Difference:    0.00000828 
Iteration:   22  Stress (raw):   0.00359331  Difference:    0.00000716 
Iteration:   23  Stress (raw):   0.00358711  Difference:    0.00000620 
Iteration:   24  Stress (raw):   0.00358173  Difference:    0.00000538 
Iteration:   25  Stress (raw):   0.00357706  Difference:    0.00000467 
Iteration:   26  Stress (raw):   0.00357299  Difference:    0.00000407 
Iteration:   27  Stress (raw):   0.00356944  Difference:    0.00000355 
Iteration:   28  Stress (raw):   0.00356633  Difference:    0.00000311 
Iteration:   29  Stress (raw):   0.00356361  Difference:    0.00000272 
Iteration:   30  Stress (raw):   0.00356122  Difference:    0.00000238 
Iteration:   31  Stress (raw):   0.00355913  Difference:    0.00000209 
Iteration:   32  Stress (raw):   0.00355730  Difference:    0.00000184 
Iteration:   33  Stress (raw):   0.00355568  Difference:    0.00000161 
Iteration:   34  Stress (raw):   0.00355427  Difference:    0.00000142 
Iteration:   35  Stress (raw):   0.00355302  Difference:    0.00000125 
Iteration:   36  Stress (raw):   0.00355193  Difference:    0.00000110 
Iteration:   37  Stress (raw):   0.00355096  Difference:    0.00000096 

plot(MDSmetrique, sub=paste("Stress, ", round(MDSmetrique$stress, digits=2)))

Non métrique

#J'affiche le stress avec le paramètre verbose=T(RUE)
MDSnonMetrique = mds(dist(t(freqs_rel_mfw), method = "euclid"), ndim=2, type="interval",verbose=T)

Iteration:    1  Stress (raw):   0.02349141  Difference:    0.05037561 
Iteration:    2  Stress (raw):   0.01787966  Difference:    0.00561176 
Iteration:    3  Stress (raw):   0.01617650  Difference:    0.00170315 
Iteration:    4  Stress (raw):   0.01541740  Difference:    0.00075910 
Iteration:    5  Stress (raw):   0.01499855  Difference:    0.00041885 
Iteration:    6  Stress (raw):   0.01474414  Difference:    0.00025441 
Iteration:    7  Stress (raw):   0.01458330  Difference:    0.00016085 
Iteration:    8  Stress (raw):   0.01447937  Difference:    0.00010393 
Iteration:    9  Stress (raw):   0.01441109  Difference:    0.00006828 
Iteration:   10  Stress (raw):   0.01436553  Difference:    0.00004556 
Iteration:   11  Stress (raw):   0.01433468  Difference:    0.00003085 
Iteration:   12  Stress (raw):   0.01431351  Difference:    0.00002117 
Iteration:   13  Stress (raw):   0.01429880  Difference:    0.00001471 
Iteration:   14  Stress (raw):   0.01428848  Difference:    0.00001033 
Iteration:   15  Stress (raw):   0.01428115  Difference:    0.00000732 
Iteration:   16  Stress (raw):   0.01427593  Difference:    0.00000523 
Iteration:   17  Stress (raw):   0.01427217  Difference:    0.00000376 
Iteration:   18  Stress (raw):   0.01426945  Difference:    0.00000272 
Iteration:   19  Stress (raw):   0.01426747  Difference:    0.00000197 
Iteration:   20  Stress (raw):   0.01426603  Difference:    0.00000144 
Iteration:   21  Stress (raw):   0.01426498  Difference:    0.00000105 
Iteration:   22  Stress (raw):   0.01426421  Difference:    0.00000077 

plot(MDSnonMetrique, sub=paste("Stress, ", round(MDSmetrique$stress, digits=2)))